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Abstract 

Two exemplary exothermic processes, synthesis of nitroglycerine in a continuous stirred tank re- 
actor (CSTR) and synthesis of the explosive RDX in a CSTR, are used to demonstrate the dangers 
of ignoring the system dynamics when defining criteria for thermal criticality or runaway. Stabil- 
ity analyses are necessary to prescribe such criteria, and for these systems prove the presence of 
dangerous oscillatory thermal instability which cannot be detected using the steady state thermal 
balances. 
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1. Introduction 

This short communication is an adjuration for the correct use of mathematical stability analysis 
in specifying critical runaway conditions for open thermoreactive systems, and is presented in the 
interests of process safety. In considering the thermal stability of such a system, operating at a 
given steady state, the first and most important question we should ask is 

Question 1. Will a small perturbation to the temperature grow uncontrollably, or decay harmlessly 
back to the steady state? 

Open dissipative dynamical systems in general may become unstable at a turning point (also 
called a limit point or a saddle-node bifurcation) or at a Hopf bifurcation to an oscillatory state, 
and this is the case for open thermoreactive systems too. Some open reacting systems in which 
thermal oscillations have been observed experimentally since 1969 are listed in Table 1. Typically 
the setup is a continuous stirred tank reactor (CSTR), which is an exemplary open system for 
observing bistable and oscillatory phenomena as well as being widely used in industry. 
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Table 1: Experimental CSTR systems in which thermal oscillations have been observed. 



Reaction 

Vapour-phase chlorination of methyl chloride 

Hydrolysis of acetyl chloride in acetone-water solvent 

Sodium thiosulfate and hydrogen peroxide reaction in aqueous solution 

Decomposition of hydrogen peroxide by Fe^"^ in HNO3 

Acid catalysed hydration of 2,3-epoxypropanol-l to glycerol 

H2SO4 catalysed hydrolysis of acetic anhydride 

Hydrolysis of methyl isocyanate (modelled, using experimental data) 



Reference 
Bush (1969) 
Baccaroetal. (1970) 
Chang and Schmitz (1975) 

Wirges (1980) 

Heemskerk et al. (1980); Ver- 
meulen and Fortuin (1986a,b) 
Haldar and Phaneswara Rao 
(1991); Jayakumar et al. (2011) 
Ball (2011) 



Oscillatory instabilities are also well-understood mathematically. The mathematical frame- 
work for identifying and characterising such instabilities in dynamical systems was set out by 
Maxwell (1868) and is presented in many excellent modern texts; rather than cite them here or 
reproduce the mathematical details the reader is referred to the expository and resources chapter 
by Ball and Holmes (2007). This mathematics of dynamical systems has long been accessible to 
chemical engineers. It has been applied to reactive thermal systems by many authors and their 
results published widely in the chemical engineering literature since the late 1950s. Some notable 
examples are as follows: Aris and Amundson (1958); Bush (1969); Gray (1969), where stability 
conditions are worked through in detail, and the abstract states explicitly "The critical condition 
for an oscillatory stable steady state is of an entirely new type (in explosion theory) not reducible 
to a tangency condition"; and Gray (1975). 

The alarm bells ring, therefore, when we see recent works specifying thermal runaway criteria 
where the the process dynamics and mathematical stability have been completely disregarded. 
Given the potential of oscillatory instability to cause serious thermal hazard, initiating thermal 
runaway, unintended thermal explosions, or worse, we feel that it is a matter of urgency to bring 
the importance of stability analysis to the attention of those sectors of the chemical engineering 
community which deal with reactive thermal hazards and runaway criteria, and correct certain 
misconceptions prevalent in the current literature regarding thermal stability criteria. This is best 
achieved by use of concrete examples of real systems. The following two published examples 
(Lu et al., 2008, 2005) were chosen as cases that graphically illustrate the pitfalls of ignoring 
CSTR dynamics, and where stability analysis reveals oscillatory thermal instability of dangerous 
amplitude that cannot be detected using classical (Semenov) ignition theory. 



2. Example I: Synthesis of nitroglycerine in a CSTR 

Lu et al. (2008) developed CSTR operating criteria for the esterification of glycerol with nitric 
acid to produce nitroglycerine, 

H2SO4 

C3H5(OH)3 + 3HNO3 C3H5(ON02)3 + 3H2O, (1) 

that were claimed to be safe, based on applying classical ignition criteria to a steady state thermal 
balance. (But see the derivation of CSTR equations given in Appendix B.) Their mass and 
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enthalpy balances for this system, expressed in dynamical form, are 

V— = VAe^/^^(CG,f - x)"(CN,f - ^XT + i^(CG,0 - CG,f - x) (2) 

at 

VCvoi— = y ((-A//°) + ^Cp {T - r°)) Ae^/^^(cG,f - xficns - ^xT - C,,iF(T - Tt) - L(T - T,). 

(3) 

Notation and quantities are defined in Appendix A. (Lu et al. (2008) took cg,o = Cgs and Tf = T° .) 

Equations (2) and (3) are not amenable to exact analysis, but the stability can be mapped by 
numerical computation. One such map is rendered in Fig. 1, where the loci of critical points are 
projected on the Ta-F parameter plane, and the values of the other system parameters are those 
given in Lu et al. (2008). To obtain the data for the curves in the figure we computed the steady 
states of Eqs (2) and (3) for increments in T^, beginning at = 230 K and with AT^ = 0.04 K, and 
for F = 141/min. At each steady state the corresponding linear eigenvalue problem was solved 
and the points where an eigenvalue changed sign were flagged. The system becomes unstable if 
such a sign change is positive. These flagged points are the bifurcation points: if a real eigenvalue 
becomes positive there is a saddle-node bifurcation, or classical ignition point, and if the real parts 
of a complex conjugate pair become positive the system exhibits a Hopf bifurcation and the onset 
of limit cycle oscillations. Then, F was allowed to vary as well as and the loci of the Hopf and 
saddle-node bifurcations was computed over Ta-F. 
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Figure 1: Loci of Hopf (solid line, magenta) and saddle-node (dashed line, blue) bifurcations in the system (2) and (3). 
The X marks a steady state considered by Lu et al. (2008). L = 80079 W/K, Tf = 273 K, Cvoi = 2703 J K ^L', AH° = 
lOOkJ/mol, ACp = -25JK-imori, n = 0.935, m = 1.117, A = 9.78 x 1022(l/mol)i-"-'«min-i, E = 122kJ/mol, 
CG.o = CG,f = 2.4904 mol/1, cn,o = 8.5535 mol/1. 

In Fig. 1 the solid (magenta) line is the locus of the Hopf bifurcations. Within the region that is 
almost enclosed by this line, and also near to the upper section of the line, the system will develop 
self-sustained thermal oscillations, regardless of the initial conditions. In the region enclosed by 
the dashed (blue) line the system has two stable steady states, between which is an unstable steady 
state. 
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Figure 2: Computed time series for selected values of Ta and F within the unstable regime in Fig. 1. (a) = 250 K, 
F - 5 1/min. (b) = 246 K, F - 9 1/min. Other parameters as given for Fig. 1 . 

The X in Fig. 1 marks one of the steady states that was considered by Lu et al. (2008). As 
can be seen, it is outside the region of oscillatory instability and the region of multiplicity and a 
check of the computed data finds the steady state system temperature T = 281 K. But if Tg, the 
coolant temperature, is reduced below 257 K (with F fixed) the Hopf bifurcation locus is crossed 
and thermal oscillations set in. A time series computed at = 250 K is shown in Fig. 2(a); the 
temperature amplitude maximum of ~295 K may cause some concern. This amplitude actually 
increases with increasing F and decreasing the time series shown in Fig. 2(b) has a seriously 
dangerous amplitude maximum. 

We would like to emphasize that the presence of these oscillatory states can only be detected 
by applying stability analysis correctly to the full dynamical system, Eqs (2) and (3). The steady 
state thermal balance as used by Lu et al. (2008) is blind to oscillatory solutions because it cannot 
answer, or even ask. Question 1 . At best a steady state thermal balance can give the curve of saddle- 
node bifurcations in Fig. 1, but since it cannot make any statement about the system stability at 
such a point it is, strictly speaking, inadequate even for specifying classical ignition/extinction 
criteria. Furthermore, a classical ignition point does not automatically infer thermal runaway or 
criticality. A check of the computed data for F = 14 1/min finds the steady state temperature 
T = 267 K at the putative 'ignition' point and T = 257 K at the 'extinction' point. These are very 
low temperatures for this system, and in this case the saddle-node locus is not useful for defining 
thermal criticality. 

We also point out that the presence of oscillatory instability may have serious implications for 
start-up and shut-down procedures. Figure 3 shows the time series computed for a slow drift or 
tuning of T^, = +0.07 K/min. In both cases the onset of thermal oscillations is violent and the 
amplitude is dangerously high. 

3. Example II: Synthesis of RDX in a CSTR 

Synthesis of the explosive cyclo-l,3,5-trimethylene-2,4,6-trinitramine, commonly known as 
RDX or hexogen or cyclonite, by reaction of hexamine with excess nitric acid is exothermic and 
known to exhibit thermal instability (Luo et al., 2002). The synthesis is often carried out in 100% 
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Figure 3: Computed time series for F - 81/min with simulated slow drift or tuning of ambient temperature, (a) 
- 0.07 K/min, (b) Ta = -0.07 K/min. Other parameters as given for Fig. 1. 



HNO3 at -20°C (Singh et al., 1999). The following overall reaction includes a secondary reaction: 
2C6H12N4 + 26HNO3 ^ CsHeNgOfi + 9CH2(ON02)2 + 5NH4NO3 + 3H2O. (4) 

RDX 

The primary reaction is 

C6H12N4 + IOHNO3 ^ C3H6N6O6 + 3CH2(ON02)2 + NH4NO3 + 3H2O, (5) 

RDX 

This nitration was modelled by Lu et al. (2005) as a CSTR process using reaction (4). Taking 
the steady state enthalpy balance they applied Semenov theory to map the classical ignition and 
extinction points over the parameter space. However, consideration of the stability of the steady 
states yields a different picture. The CSTR mass and enthalpy balances used by Lu et al. (2005) 
(but see Appendix B), expressed in dynamical form, are 

-VAe'^'^^c^c" + F{cf - c) (6) 
V (i-AH°) + ACp (T - D) Ae^/^^CNC" - CvoiF(r - Tf) - L{T - T,). (7) 

Notation and quantities are defined in Appendix A. 

We carried out a numerical stability analysis of this system using a procedure exactly analogous 
to that described in section 2, and found that it, too, exhibits oscillatory instability. In Fig. 4 the 
solid (magenta) line marks the locus of the Hopf bifurcations dashed (blue) line marks the locus 
of saddle-node bifurcations, or classical ignition points F-T^, for the given values of the other 
parameters. The qualitative similarity to Fig. 1 may be noted, and the interpretation given in 
section 2 may be applied here too. 

The value of Ta and of F marked by the X in Fig. 4 were chosen to compute a time series; this 
is shown in Fig. 5. The high amplitude of the thermal oscillation, at 317 K, is of concern: above 
~300K exothermic side-reactions are reported to take over (Hale, 1925), rendering the system 
extremely unsafe. 
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Figure 4: Loci of Hopf bifurcations (solid line) and saddle-node bifurcations in the system (6) and (7). X marks the 
values of T.^ and F for which the time series in Fig. 5 was computed. L - 50 W/K, Tf - 273 K, Cyoi - 2940 JK"'l"\ 
AH = 90.938 kJ/mol, ^Cp = -322 JK-'mor', n = 1.28, A = 1.234 x 10^(l/mol)"min-i, E = 47.152kJ/mol, 
Cf = 4.805 mol/1, cn =21.35 mol/1, T° = 298 K. 




Figure 5: A time series integrated from Eqs (6) and (7), with T-^ - 254 K, F — 0.0185 1/min and other parameters as 
given for Fig. 4. 



The standard reaction enthalpy includes the enthalpy of mixing nitric acid and hexamine. From 
standard enthalpies of formation Luo et al. (2002) found A//°(298 K)= -90.938 kJ/mol hexamine 
for the overall reaction (4), and this is the value we used in the computations. However, for the 
main reaction (5) A//°(298 K) = - 153.31 8 kJ/mol hexamine for the reaction only, from standard 
enthalpies of formation. Both of these values are at odds with that measured by Dunning et al. 
(1952) using straightforward calorimetry, at least 64kcal/mol= 270kJ/mol hexamine at -35.5°C 
in 99% nitric acid, which gives A//°(298 K)= 290kJ/mol hexamine; presumably exothermic sec- 
ondary reactions were involved. In any case, we can expect the oscillatory instabilities to be more 
violent and to occur over more of the parameter space for a more exothermic process. 

It is helpful to study the computed solutions rendered in terms of the temperature T over in 
Fig. 6, for selected values of F. Stable steady states are marked with solid lines (blue) and unstable 
steady states are marked with dashed lines (red). Also marked are the maximum and minimum 
amplitude envelopes of the limit cycles (which were also computed and analysed for stability), 
with thick dotted lines (magenta), from their origin at a Hopf bifurcation to their termination at 
a second Hopf bifurcation (in (c) and (d)) or at a homoclinic bifurcation (in (a) and (b)). These 
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Figure 6: Stable (solid lines, blue) and unstable (dashed lines, red) steady states and amplitude envelopes of limit 
cycles (thick dotted lines, magenta) in the system (6) and (7). (a) F - 0.025 1/min, (b) F - 0.020 1/min, (c) F - 
0.0185 1/min, (d) F - 0.017 1/min. Other parameters as given for Fig. 4. 



bifurcation diagrams are a rich source of information about the thermal stability of the system. 

In Fig. 6(a) we see that classical ignition occurs at the lower saddle-node bifurcation to a 
stable steady state on the upper stable solution branch. This does not necessarily represent thermal 
runaway, although the temperature of ~313K may be too high for comfort. The temperature is 
reduced by quasistatic reduction of — but as the Hopf bifurcation is crossed rapidly growing 
thermal oscillations set in, and at the terminus T ~ 322 K. We also see that extinction does not 
occur at the second saddle-node bifurcation, as assumed by Lu et al. (2005), but at this oscillatory 
terminus (provided thermal runaway has not already set in). 

In Fig. 6(b) the ignition is non-classsical; it occurs from the lower saddle-node bifurcation to 
a limit cycle that grows in amplitude with Ta, reaches a maximum, then declines. 

Steady state multiplicity is barely present at the value of F for which Fig. 6(c) was computed. 
Non-classical ignition from the lower saddle-node bifurcation occurs to a limit cycle, but extinction 
in this case is classical, occuring at the upper saddle-node bifurcation. 

The system globally has no multiplicity in Fig. 6(d) and thermal runaway, if it occurs, must be 
the result of oscillatory behaviour. 



4. Discussion and summary 

The two examples treated here provide ample illustration of the dangers inherent in using a 
steady state thermal balance and ignoring the dynamics when determining thermally safe operating 
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criteria for open thermoreactive processes. Classical (Semenov) ignition theory was originally 
applied to explosive solid materials then extended to reacting gases at low pressure, for which it 
works well. No reference is made in Semenov theory to oscillatory solutions, because reactant 
depletion is ignored (the infinite pool approximation) and periodic solutions are automatically 
forbidden in a system evolving in just one dynamical variable — the temperature. Semenov theory 
is insufficient for analysing open flow systems with reactant depletion, for which the steady states 
and their stability must be fully characterised. For the simple case of a single first-order (or 
pseudo first-order) exothermic conversion in a CSTR the global point of onset of multiplicity has 
been calculated analytically (Ball, 1999) and analysis combined with numerics has characterised 
the oscillatory states and identified thermal runaway due to the hard onset of oscillatory instability 
at a subcritical Hopf bifurcation (Ball and Gray, 1995; Ball, 201 1). For more complicated or more 
accurately modelled systems numerical eigenvalue analyses can be carried out, as we have done 
here. 

More broadly, we are concerned that crucial knowledge that was available (and, of course, is 
still available) to and made good use of by chemical engineers from the 1950s through the 1970s 
seems to have been forgotten or never was learned in some relevant quarters. StabiHty theory and 
mathematics in the context of thermoreactive systems has nothing to do with steady state thermal 
balances and everything to do with protecting human Uves by assaying thermal runaway conditions 
correctly and accurately. 

We summarize the message of this communication as follows: 

• It is incorrect, dangerous, and inconsistent with the existing literature on this subject to 
assume that thermal runaway is governed only by classical ignition, and use a steady state 
enthalpy balance to answer the wrong question: What are the boundaries of classical ignition 
and extinction points over the parameter space? 

• The current study has shown that 

- classical ignition need not be accompanied by thermal runaway, 

- thermal runaway may be non-classical, with ignition to an oscillatory state, 

- non-classical thermal extinction may occur from an oscillatory state, 

- oscillatory thermal runaway may occur in the complete absence of steady state multi- 
plicity, and 

- oscillatory thermal instability is endemic to these systems; it is real, it is not some 
numerical or other artefact, it may occur with violent abruptness, it is potentially dan- 
gerous. It must be taken into account in devising thermal runaway or safe operational 
criteria. 

• The right question to ask of open thermoreactive systems is Question 1 : Given a steady state, 
will a small perturbation grow or decay? Answering this question necessarily involves the 
dynamics; i.e., a stability analysis using the well-founded mathematics of stability theory. 

In view of these results, we strongly recommend that numerical stability analyses be carried out, 
over the relevant system and operating parameter space, for all exothermic reactions that are run 
in a CSTR, using procedures that have been in the pubHshed literature for many decades. 
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Appendix A. 



Table A.2: Notation and definitions. 



Symbol or 


Definition 


quantity 




A 


pre-exponential factor 


c 


concentration of reactant (glycerol or hexamine, mol/1) 




molar heat capacity (J/(mol K) 


C'vol 


average volumetric specific heat (J T^K"') 


ACp 


reaction heat capacity (J/(mol K) 


E 


activation energy (J/mol) 


F 


volumetric inflow rate ( 1/min) 


H 


enthalpy of formation (J/mol) 


AH 


reaction enthalpy (J/mol) 


k 


kiT), temperature-dependent rate constant 


L 


combined heat transfer coefficient (W/K) 


n, m 


reaction orders 


t 


time (min) 


T 


temperature (K) 


V 


volume of reaction mixture (1) 


X 


fractional conversion of reactant 


Subscripts, superscripts, overscripts 





initial value of quantity 


o 


standard state 




time-changing quantity 


a 


of the coolant 


A,B, S 


of reactant A, B, of solvent S 


f 


of the feed stream 


G 


of glycerol 


N 


of nitric acid 



Appendix B. 

The enthalpy balance appears to have been given incorrectly in Lu et al. (2005) and Lu et al. 
(2008), with the reaction enthalpy AH° being counted twice. The correct enthalpy balance for the 
adiabatic CSTR housing a first order exothermic reaction A ^ B in a solvent S was given in Ball 
(1999); it simply expresses the law of energy conservation and is reproduced here for convenient 
reference: 

[caHa(T) + (Cf - Ca) Hj,(T) + csHsiT)] = 

at 

F [cfHAiTf) - caHa(T) - (Cf - c)Hb(T)] 

+ F[cs,fHs(Tt)-csHs(T)] (B.l) 

Using the thermodynamic relation Cp = (dH/dT)p Eq. 6 may be developed into 

[CT + (-AH) ca] = -F [{CT + (-AH) ca) - {CTf + (-AH) Cf)] . (B.2) 
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The differential mass balance 

= -Vk(T)c^ + F(cf - ca) (B.3) 

at 

is substituted into Eq. B.2 to give, for arbitrary initial conditions, 

= V(-AH)k(T)cA - FC(T - Tf). (B.4) 

at 

One may include any number of reaction components and heat loss terms in the fundamental 
energy conservation expression Eq B.l. The correct enthalpy balance, Eq. B.4 emerges naturally 
by considering the dynamics; i.e., the rate of enthalpy generation, the left hand side of Eq. B.l. 
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